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ABSTRACT 


Computer  simulation  of  the  mechanics  of  molecular  systems 
is  a  popular  and  powerful  method  for  understanding  chemical 
processes  The  complexity  of  modeled  chemical  systems  has 
advanced  from  hard  spheres  and  rare  gases  to  liquid  solutions  and 
hiomolecules  Such  simulations  are  computationally  intensive  and 
thus  are  limited  by  the  speed  of  available  computers  This  paper 
describes  the  use  of  specialized  hardware,  a  high  speed  floating 
point  array  processor  (AP),  to  dramatically  speed  up  molecular 
mechanics,  in  other  words  molecular  dynamics.  Monte  Carlo,  and 
energy  minimization  calculations  Although  the  array  processor  is 
a  cost  effective  solution  for  computationally  intensive  problems  in 
terms  of  hardware  (full-time  AP  usage  is  equivalent  to  two  to  eight 
hours  per  day  of  Cray-1  time),  Its  full  speed  comes  at  the  expense 
of  programming  in  a  relatively  difficult  parallel  assembly  language. 
Since  the  architecture  of  the  machine  is  dramatically  different 
from  conventional  computers  and  utilizing  its  fast  speed  necessi¬ 
tates  using  this  architecture  on  the  assembly  language  level,  the 
proper  design  and  implementation  of  algorithms  is  critical  The 
molecular  mechanics  software  design  discussed  here,  consisting  or 
12  000  lines  of  C  and  7  000  lines  of  AP  assembly  language  code,  is 
quite  genera]  and  has  been  used  to  study  systems  ranging  from 
rare  gases  to  biomolecules  This  implementation  yields  effective 
speeds  approximately  thirty-five  times  faster  than  a  dedicated  DEC 
VAX  n/780  computer  with  floating  point  accelerator  and  optimized 
VMS  Fortran,  thus  allowing  simulations  to  be  run  in  a  week  and  a 
half  on  the  AP  which  would  require  a  year  of  dedicated  VAX  time. 
The  flexibility  of  the  IDOC  operating  system,  whose  source  code  is 
accessible  and  can  be  modified  to  optimize  performance,  com¬ 
bined  with  the  modern  features  of  the  C  language,  have  made  this 
Implementation  much  easier,  by  providing  a  convenient  and  power¬ 
ful  environment  in  which  to  imbed  the  hand-coded  AP  assembly 
language  modules.  Applications  to  date  range  from  the  molecular  j 
dynamic  calculation  of  infrared,  Raman,  and  electronic  spectra  in  ; 
gas  and  liquid  solutions  to  the  calculation  of  thermodynamic  quan-  j 
tJties  for  water  and  the  simulation  of  the  molecular  dynamics  of 
solution  reactions  and  of  polypeptides. 
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I.  INTRODUCTION 

-In  recent  years  molecular  mechanics,  the  computer  simulation  of  molecu¬ 
lar  systems  using  molecular  dynamics,  Monte  Carlo,  and  energy  minimization, 
has  eme-ged  as  a  powerful  tool  for  investigating  and  understanding  chemical 
properties  and  processes,.  While  straightforward  in  principle,  these  techniques 
can  be  unusually  demanding  computationally  and  thus  the  size  and  complexity 
of  systems  which  can  feasibly  be  simulated  is  determined  by  the  speed  and  avai¬ 
lability  of  computer  hardware /?ln  this  paper  we  discuss  a  particular  solution  to 
the  computational  needs  of  molecular  mechanics,  the  use  of  specialized 
hardware,  a  high  speed  array  processor,  in  our  case  a  Floating  Point  Systems, 
Inc  AP-i209f ’whose  use  we  pioneered  after  purchasing  serial  number  two  /  In 
other  papersh^we  have  also  discussed  an  alternative  solution,  the  division  of  the 
problem  among  an  array  of  different  processors  operating  in  parallel  Yet  a 
third  solution*is  to  use  a  vector  processing  machine  such  as  a  Cray-1  We  will 
focus  here  primarily  on  molecular  dynamics  within  the  array  processor  program 
package  for  molecular  mechanics  we  have,  developed,  called  ^Newton*  whose 
OOoception  is  described  in  earlier  papers.1-*^!  section  II  we  examine  the  archi¬ 
tecture  of  the  AP-i20BfJas  this  is  necessary  to  understand  how  to  efficiently  use 
the  machine  >»ln  section  III  we  lay  out  the  specifics  of  the  structure  of  the  pro¬ 
gram  package  we  have  developed  for  molecular  dynamics  Finally,  in  section  IV 
we  present  the  results,  analyze  potential  array  processor  improvements,  and 

Cint  out  the  advantages  and  importance  of  the  environment  in  terms  of 
iguage  and  operating  system* 

Classical  molecular  dynamj:s*  is  the  simulation  of  systems  by  the  numeri¬ 
cal  integration  of.  for  example,  Newton’s  Second  Law,  to  obtain  particle  trajec¬ 
tories,  l.e  coordinates  and  momenta  as  functions  of  time,  and  then  the  compu¬ 
tation  of  physical  properties  as  averages  over  these  trajectories  In  the  Monte 
Carlo  technique,  particles  are  moved  artificially  as  the  result  of  random  number 
generation,  rather  than  dynamically.  The  Monte  Carlo  technique  is  less  general 
than  molecular  dynamics  in  that  it  is  appropriate  for  the  calculation  of  proper¬ 
ties  which  depend  upon  averages  over  coordinates,  but  not  for  properties  which 
(are  explicitly  time  dependent  These  fields  are  reviewed  by  Alder,5  Wood  and 
Erpenbeck.8  Valleau  and  Whittington,7  McDonald,8  Ifinder,®  and  Wood10  and  are 
the  subject  of  a  workshop,11  Energy  minimization12-15  is  the  search  for  the 
minimum  energy  of  the  system  as  a  function  of  particle  coordinates. 

Early  liquid  state  molecular  dynamic  simulation18  Involved  hard  sphere 
atoms.  Later,  soft  potentials,  eg  Lennard-Jones.  were  introduced.17- 18 
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Computer  simulation  became  a  too!  for  testing  approximate  theones  for  atomic 
liquids  16  More  recently,  with  increasingly  powerful  software  and  faster 
hardware,  attention  has  focused  on  more  complex  molecular  systems,  ranging 
from  water’0  to  proteins20  and  nucleic  acids 

Array  processors  can  bring  much  greater  computer  power  to  bear  on  chem¬ 
ical  problems  than  was  previously  available  2 21  ■*2  However,  their  use  carries 
the  serious  disadvantage  that,  at  least  up  to  the  present,  their  maximum 
efficiency  has  only  been  gained  by  assembly  language  programming  as  an 
efficient  compiler  for  such  unusual  architecture  is  not  yet  available  Thus  the 
development  of  efficient  and  general  purpose  software  for  array  processors  is 
very  important  as,  in  the  long  run  software  costs  will  usually  dominate 
hardware  costs  even  in  a  University  environment.  We  thus  hope  that  the  tech¬ 
niques  discussed  here  which  we  have  used  in  the  development  of  Newton,  our 
tool  for  molecular  mechanics  on  an  array  processor,  will  prove  useful  to  others 
using  array  processors  for  molecular  mechanics  as  well  as  for  other  classes  of 
problems. 

D.  ARRAY  PROCESSOR 

There  are  three  ways  to  build  faster  computer  hardware.  The  first  is  to 
Increase  the  speed  of  the  logical  elements  themselves,  the  second  is  to  horizon¬ 
tally  spread  out  the  calculation  among  many  parallel  elements  all  working  at 
once,  and  the  third  is  to  vertically  spread  out  the  steps  of  the  calculation  in 
time  along  a  pipeline,  so  that  many  successive  different  calculations  can  trickle 
down  the  pipeline  together  Arrays  of  processors1' 2  are  an  example  of  highly 
parallel  architecture  and  supercomputers  such  as  the  Cray-1  and  Cyber  205 
make  extensive  use  of  pipelines  (as  well  as  parallelism). 

The  Floating  Point  Systems  AP- 1203  (AP)  array  processor  is  a  special  purpose 
computer,  which  combines  all  three  approaches,  logic  speed,  parallelism,  and 
pipelining  but  pushes  none  of  these  to  its  extreme  limit.  It  is  designed  for  very 
fast  processing  of  floating  point  arithmetic  23  In  comparison,  logical  instructions 
are  much  slower  and  more  cumbersome.  It  is  not  an  array  of  processors  as 
many  infer  from  the  name,  but  a  single  parallel  and  pipelined  processor  whose 
architecture  is  easily  exploited  for  array  type  operations  The  maximum  floating 
point  rate  is  12  million  operations  per  second  (12  megaflops)  This  speed  puts 
the  AP  in  the  same  computational  class  as  many  larger  main  frame  computers, 
as  indicated  in  Table  1  Because  it  is  a  peripheral  processor,  the  AP  requires  a 
host  computer  to  initiate  data  and  prog  ram  transfers  to  and  from  it  The  nature 
of  the  host  interface  as  well  as  the  host  operating  system  are  thus  important 
factors,  as  will  be  discussed  below 

Why  Use  an  Array  Processor? 

The  major  reason  that  the  use  of  an  array  processor  is  attractive  is  its  cost 
effectiveness  This  is  shown  in  Table  I  from  the  work  of  Bucy  and  Senne24  as  is 
described  in  detail  by  Karplus  and  Cohen.25  Columns  one  and  two  show  comput¬ 
ers  that  have  been  commonly  used  in  scientific  numerical  applications  and  their 
theoretical  speed,  as  measured  in  millions  of  floating  point  operations  per 
•econd  (megaflops)  Columns  three  and  four  show  the  performance  in  terms  of 
time  and  achieved  megaflops  for  a  sample  application  Columns  ftve  and  6ix 
•how  the  motivating  economic  reason  for  using  the  ap-izcb  It  has  the  lowest  cost 

Sr  achieved  megeflop  and  one  of  the  lowest  overall  installation  costs  However 
e  last  column  shows  where  the  real  price  is  paid.  The  AP- 1»3  has  a  substan¬ 
tially  higher  program  development  time,  in  addition  to  a  costly  learning  curve  as 
users  must  become  familiar  with  its  unusual  and  challenging  parallel  assembly 
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Programming  Techniques 

Three  methods  of  programming  the  array  processor  are  available  The  level 
highest  in  abstraction,  but  slowest  in  speed,  is  the  Fortran  compiler.  Fortran  is 
a  language  familiar  to  most  computational  chemists,  but  unfortunately  no  exist¬ 
ing  Fortran  compiler®6"28  for  the  AP-1203  generates  efficient  or  compact  code. 
Compiled  code  is  thus  lengthy  (one  can  quickly  consume  more  program  source 
memory  than  is  available  for  the  AP-1203)  and  therefore  necessarily  slow,  since 
the  machine  is  synchronous  Most  applications  would  thus  require  the  storage  of 
instructions  in  main  data  memory  which  is  costly  not  only  due  to  the  loss  of 
main  data  memory  (two  locations  per  each  instruction  stored)  but  also  in  the 
overhead  of  the  overlaying  process  both  in  terms  of  the  bookkeeping  required 
and  the  actual  transfer  of  the  instructions  to  and  from  main  data  memory  as  the 
overlaying  is  done.  The  improvement  in  program  speed  and  performance  of  For¬ 
tran  in  the  AP-1203  over  a  VAX  n/7B0  can  be  as  little  as  a  factor  of  two  or  three  21 
This  degrades  even  further  if  many  data  transfers  to  and  from  the  host  or  much 
loading  of  program  source  memory  from  main  data  memory  is  required 

The  second  level  of  usage  is  to  vectorize  the  algorithm  into  a  form  in  which 
the  problem  can  be  solved  by  a  series  of  calls  to  canned  vector  routines  written 
In  assembly  language  by  floating  Point  Systems  (FPS)  This  can  be  an  improve¬ 
ment  over  the  Fortran  compiler,  but  since  the  code  for  the  individual  vector 
operations  cannot  be  overlapped,  the  full  advantage  cannot  be  taken  of  the 
speed  of  the  AP.  Two  modes  of  operation  are  possible,  one  in  which  the  canned 
routines  are  called  on  an  individual  basis  from  the  host  C  or  Fortran  program 
and  the  other  where  these  cedis  are  linked  together  in  a  rudimentary  higher 
level  language,  know  as  the  vector  function  chainer.  The  vector  function  pro¬ 
gram  is  then  compiled  into  AP  assembly  code,  and  executes  the  set  of  opera¬ 
tions  as  a  group.  It  Is  vitally  important  to  use  the  vector  function  chainer  to 
string  these  operations  together  so  that  the  overhead  involved  in  loading  and 
starting  the  AP  is  only  paid  once  However,  the  vector  function  chainer  gen¬ 
erates  even  worse  code  than  the  Fortran  compiler,  so  any  loops  coded  in  the 
vector  function  program  can  also  degrade  performance  A  typical  Increase  in 
speed  over  that  of  a  VAX  ti/7B0  is  three  to  twenty  using  this  approach,  as  it  is  very 
dependent  on  how  well  the  particular  canned  functions  encompass  the  problem 
being  solved  Molecular  mechanics  calculations  tend  to  be  at  the  low  end  of  this 
range 

The  third  level  of  usage,  needed  to  realize  the  full  potential  of  the  AP,  is  to 
hand-code  commonly  used  algorithms  In  AJM20B  assembly  language  This 
presents  a  problem  to  most  computational  chemists  in  that  it  may  require  skills 
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unfamiliar  to  them  In  order  to  investigate  the  reason  for  the  necessity  and 
Inherent  difficulty  of  assembly  language  coding,  a  knowledge  of  the  architecture 
of  the  AP  is  helpful 

Architecture 

The  array  processor  is  a  parallel  and  pipelined  machine  26  A  parallel  com¬ 
puter  is  one  that  can  perform  more  than  one  operation  in  a  single  instruction 
The  array  processor  consists  of  several  independent  memories  and  arithmetic 
units,  as  ehown  in  Pig  1.  all  of  which  can  be  addressed  or  initiated  within  a  sin¬ 
gle  instruction  of  the  AP-1203  There  are  three  types  of  memory.  Program 
source  memory,  which  is  64  bits  wide,  is  used  to  store  assembly  language 
instructions  Data  is  not  stored  intermixed  with  instructions,  as  in  most  com¬ 
puters.  but  resides  primarily  in  main  data  memory  (36  bits  wide)  which  is  avail¬ 
able  in  two  speeds,  known  as  "feist"  and  "slow"  memory  Data  can  also  be  stored 
in  table  memory  which  can  consist  of  both  Head  Only  Memory  (ROM)  and  Random 
Access  Memory  (RAM)  In  our  configuration  we  have  64  K  words  of  main  data 
memory.  2K  words  of  table  ROM,  4  K  words  of  writable  table  memory,  and  ]  .5  K 
words  of  program  source  memory  The  actual  hardware  of  the  machine  can  sup¬ 
port  up  to  4  K  of  program  source,  4  K  of  writable  table  memory,  and  512  K  of 
main  data  memory  (however,  only  one  individual  64  K  bank  can  be  accessed  at 
any  one  time).  In  addition  there  are  two  banks,  DP^X  and  DPY,  of  thirty  two  (38  bit) 
registers  called  data  pads  (only  eight  of  which  can  be  accessed  at  a  time)  for  the 
storage  of  intermediate  results  Addresses  (pointers)  and  integer  parameters 
are  stored  in  sixteen  (16  bit)  integer  registers,  known  as  S  pads  In  addition  to 
the  memories  there  are  two  separate  floating  point  arithmetic  processors,  an 
adder  and  a  multiplier.  The  floating  adder  is  capable  of  a  variety  of  operations 
such  as  addition,  subtraction,  logical  and  as  well  as  logical  or  There  is  also  an 
Integer  arithmetic  unit  which  can  be  initiated  every  instruction  and  whose  result 
is  available  for  use  in  the  instruction  in  which  it  was  initiated 

F*  2  shows  the  pipeline®  nature  of  the  AP-1203  The  floating  multiplier  dI 
the  AP  is  a  three  stage  pipeline,  and  thus  a  multiply  operation  started  in  one 
instruction  will  take  a  minimum  of  three  instructions  to  complete  A  new  multi¬ 
ply,  however,  can  be  initiated  at  every  instruction  and  each  initiation  causes  the 
preceding  operations  to  be  pushed  through  the  stages  of  the  pipeline  If  the 
preceding  operations  are  not  pushed,  the  results  stay  in  the  pipeline  until 
pushed  out  later  by  initiating  other  multiplies.  The  floating  adder  is  a  two  stage 
pipeline  with  similar  properties  to  the  multiplier  The  memory  units  are  also 
pipelined  in  the  sense  that  a  memory  fetch  initiated  in  one  instruction  does  not 
complete  and  cannot  be  used  until  three  (or  two,  in  the  case  of  table  memory) 
instructions  later. 

Program  branching  can  only  be  performed,  with  a  few  exceptions,  on  the 
output  of  either  the  floating  adder  or  the  integer  arithmetic  unit  For  example, 
to  branch  If  a  particular  floating  number  is  negative,  the  number  must  first  be 
generated  in  the  adder,  which  requires  two  instructions  On  the  third  instruction 
the  result  is  available  from  the  adder  and  on  the  fourth,  and  not  before,  the 
result  may  be  tested  and  branched  upon.  This  is  the  principle  reason  that  cal¬ 
culations  which  require  a  significant  amount  of  logic  do  not  perform  as  well  on 
the  AP-iSOB  Branching  based  on  integer  operations  is  not  as  Blow  because  the 
result  Is  available  in  the  instruction  in  which  the  operation  is  initiated  and  may 
thus  be  tested  in  the  following  Instruction  In  either  case,  building  logic 
efficiently  into  an  AP-1209  assembly  language  program  can  be  quite  cumbersome 

Another  important  feature  of  the  architecture  of  the  AP-1203  is  that  each 
separate  memory  or  register  unit  can  have  one  read  or  write  (in  the  case  of 


PROGRAM  SOURCE  I  OPX  OPY  S  PAD 


Figure  1.  AP-1200  multiple  functional  units  Illustrating  the  parallel  nature  of  the 
array  processor  (AP)  It  consists  of  multiple  memory  and  arithmetic  units  which 
can  all  be  accessed  or  initiated  in  a  single  instruction  cycle 
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registers  like  DPX  or  DPY  both  a  read  and  a  write)  initiated  per  instruction  Thus 
the  placement  of  data  among  the  various  memories  is  extremely  important  in 
that  the  more  spread  out  the  data  the  more  data  can  be  accessed  in  an  indivi¬ 
dual  instruction  For  this  reason  it  is  important  to  have  writable  table  memory 
in  addition  to  main  data  memory  The  difficulty  encountered  with  the  use  of 
writable  table  memory  is  that  changing  data  stored  in  that  memory  is  a  multiple 
instruction  process  and  very  slow  in  comparison  with  writing  into  main  data 
memory  Whether  it  is  important  to  have  the  fast  rather  than  the  slow  main 
data  memory  is  highly  dependent  on  the  algorithm  being  coded  For  molecular 
mechanics  code,  which  is  usually  computation  rather  than  memory  bound,  the 
slow  memory  seems  adequate 

As  can  be  seen  in  Fig  3,  all  these  units  interconnect  through  a  complicated 
data  pad  bus  structure  It  is  contention  for  the  data  pad  bus,  which  can  contain 
only  one  value  per  instruction,  that  is  the  most  common  difficulty  m  A P  assem¬ 
bly  language  programming  Thus  efficient  programming  of  the  A?-)2C3  involves 
exploiting  the  parallel  nature  of  the  components  of  the  AP  and  the  inner  connec¬ 
tivity  of  the  data  paths  to  give  the  maximum  throughput  This  means  that  the 
manner  in  which  a  particular  part  of  an  algorithm  is  coded  in  the  early  part  of 
the  code  has  long  reaching  effects  in  the  later  stages  of  coding  the  algorithm  In 
many  instances  it  is  necessary  to  code  an  algorithm  iD  more  than  one  way  and 
then  examine  which  scheme  can  be  made  the  most  efficient  This  can  require  a 
high  degree  of  patience  and  persistence 

The  AP-1203  interfaces  to  the  host  computer  (in  our  case  a  DEC  tax  11/750) 
through  two  (or  possibly  three)  channels,  as  shown  in  Fig  3  The  virtual  front 
panel  is  used  by  the  host  to  control  the  AP-1203  It  does  not  exist  in  reality  but  is 
a  set  of  registers  (in  our  case  on  the  VAX  UNT3JS)  which  can  be  examined  and  set 
by  the  host  computer  The  UNDt  device  driver  for  the  AP  uses  these  registers  to 
start  and  stop  the  AP,  to  examine  and  deposit  into  memories  and  registers,  and 
to  initiate  Direct  Memory  Accesses  (DMA's)  over  the  UNIBUS  Any  data  manipula¬ 
tion  done  using  the  front  panel  is  necessarily  very  slow.  Thus,  data  and  program 
Instructions  are  normally  transferred  between  the  host  and  AP  via  the  DMA  pro¬ 
cess  which  can  occur  at  approximately  a  megabyte  per  second  Due  to  the 
difference  in  word  lengths  between  the  usual  host  floating  point  representation 
(32  bits)  and  the  AP  representation  (3B  bits),  floating  point  numbers  are  con¬ 
verted  "on-the-fly"  during  the  DMA  process  by  the  APs  interface  hardware  To 
preserve  the  full  precision  of  the  AP's  representation  requires  using  the  much 
■lower  front  panel  or  DMA'ing  out  the  data  in  successive  passes  Preserving  full 
precision  is  important  if  numbers  are  to  later  be  reloaded  back  into  the  AP  and 
the  calculation  continued,  as  the  result  will  not  be  the  same  as  would  be 
obtained  by  doing  the  same  uninterrupted  calculation  in  the  AP  if  the  numbers 
are  rounded  to  32  bits  The  AP  has  a  much  Larger  dynamic  range  (10',5S  to 
10*155)  than  the  usual  host  representation  (10-*  to  10**)  and  the  only  way  in 
which  these  larger  floating  point  numbers  can  be  extracted  from  the  AP-1203 
preserving  the  full  36  bits  of  precision  is  to  use  the  front  panel  However,  if  the 
floating  point  numbers  are  within  range  of  the  host’s  floating  point  representa¬ 
tion,  precision  can  be  preserved  by  DMA'ing  out  the  data  in  multiple  passes  30 
First  the  actual  data  is  DMA'ed  out,  undergoing  convergent  rounding  to  32  bits 
Next  the  data  is  DMA'ed  back  to  the  AP  and  subtracted  from  its  original  value  in 
the  AP  The  result  of  this  subtraction,  the  remainder.  Is  then  DMA  ed  to  the  host 
These  two  arrays  (the  rounded  data  and  remainder)  can  then  be  used  to  recon¬ 
struct  the  original  value  In  the  AP  by  DMA'ing  them  into  the  AP  separately  and 
adding  them  together  in  the  AP.  This  requires  the  allocation  of  scratch  areas  in 
the  AP's  memory. 
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Figure  3.  AP-1208  host  Interface  and  internal  connections  The  various  data  paths 
connecting  the  Individual  functional  units  Internal  to  the  AP-120H  are  shown  in 
the  right  portion  of  the  figure  On  the  left  are  shown  the  data  paths  that  inter¬ 
face  the  AH  to  the  host  computer 
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The  DMA  channel  also  is  a  bottleneck  when  large  volumes  of  data  are  to  be 
transferred  to  and  from  the  host  computer  The  AP  can  use  approximately  B0% 
of  the  bandwidth  of  the  U (Cat's  A  typical  transfer  rate  is  approximately  0  8  KB 
(megabytes  per  second)  Particularly  operating  system  dependent  is  the  sever¬ 
ity  of  overhead  paid  for  initiating  each  transfer  On  VAX's,  the  poor  hardware 
design  of  the  UK BUS  adapter  also  lowers  the  bandwidth  However,  the 
throughput  can  be  maximized  by  the  using  the  asynchronous  capability  of  the 
DUA  channel  A  third  optional  DUA  channel,  know  as  an  IOP,  operates  m  a  manner 
similar  to  the  standard  DUA  channel  except  that  it  has  less  general  format 
conversion  capabilities  and  lacks  a  sophisticated  handshaking  protocol  which 
may  make  it  inadequate  for  some  applications 

Another  important  architectural  feature  of  the  AP-120B  is  that  al!  main  data 
memory  locations  are  38  bits  and  designed  to  accommodate  a  single  floatmg 
point  number  While  it  is  a  quite  possible  (and  easy)  to  store  a  single  16  bit 
Integer  in  each  38  bit  location,  in  many  instances  it  would  be  more  efficient  to 
pack  more  than  one  integer  in  each  floating  point  word  This  would  provide  sav¬ 
ings  not  only  in  the  amount  of  memory  required  but  also  in  execution  time  by- 
cutting  down  the  number  of  memory  accesses  required  to  retrieve  the  data 
Although  the  38  bit  word,  as  shown  in  Pig.  4.  consists  of  three  different  fields,  a 
16  bit  low  mantissa,  a  12  bit  high  mantissa  and  a  10  bit  exponent,  not  all  of 
these  bits  are  accessible  to  an  AP-120S  assembly  language  program.  All  16  bits  of 
the  low  mantissa  and  all  10  bits  of  the  exponent  can  be  placed  mto  the  S  pad 
registers,  but  only  B  of  the  12  bits  of  the  high  mantissa  (the  so-called  table 
lookup  bits)  can  be  accessed  by  the  AP-1203  program  Furthermore,  there  is  no 
way  to  write  these  integers  back  mto  the  packed  word  from  an  AP  program  and 
these  packed  integers  cannot  be  loaded  by  the  DUa  channel  mto  the  AP-i203,  but 
must  be  loaded  using  the  front  panel  which  can  examine  and  deposit  all  the  bits 
of  the  low  mantissa,  high  mantissa  and  exponent  separately. 

Ill  NEWTON 

Our  tool  for  molecular  mechanics,  Newton  consists  of  a  variety  of  hardware 
processors  and  a  generalized  software  package  which  will  be  described  in  this 
section  The  software  portion  of  Newton  consists  of  12  000  lines  of  C  code  and 
7000  lines  of  AP  code  It  is  designed  to  be  a  modular  and  generalized  approach, 
applicable  to  systems  from  rare  gases  and  diatomics  to  polypeptides,  proteins 
and  nucleic  acids  in  solution  While  Newton  can  be  used  for  Vonte  Carlo  and 
energy  minimization,  Its  major  use  to  date  has  been  for  molecular  dynamics 
based  calculations,  and  it  is  these  we  will  describe  in  more  detail  The  array 
processor,  as  shown  in  Fig  3,  is  connected  to  a  VAX  u/?50  as  the  host  processor 
The  second  DMA  port  can  be  connected  to  a  PDF  n/54  which  serves  as  the  host  to  a 
three  dimensional  graphical  display  system,  an  Evans  and  Sutherland  Picture 
System  All  four  of  these  processors  can  be  activated  simultaneously  to  allow 
the  computation  of  molecular  dynamics  (and  derived  properties  such  as  spectra 
and  thermodynamics)  and  the  real  time  display  of  the  atomic  motions  or 
derived  properties  on  the  Picture  System  as  they  are  calculated 

We  know  of  four  other  applications  of  aspects  of  molecular  mechanics  which 
have  been  made  using  array  processors  Pottle,  Scheraga  and  co-workers35 
describe  a  package  for  energy  minimization  of  proteins  Although  the  problems 
of  force  evaluation  are  similar,  the  approaches  chosen  are  different  Newton  is 
not  confined  to  any  one  set  of  potential  surfaces,  such  as  the  electrostatic  plus 
Lennard-Jones  potential  surfaces  used  by  Pottle,  but  is  capable  of  implementing 
essentially  arbitrary  potential  surfaces,  as  will  be  seen  below  However,  we  have 
also  paid  the  price  for  this  generality  While  their  Inner  loop  code  operates 
using  80%  of  the  theoretical  maximum  floating  point  speed  of  the  array 
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Figure  4.  AP-1203  floating  point  representation  A  main  data  memory  word  is 
divided  into  three  parts  a  16  bit  low  mantissa  field,  a  12  bit  high  mantissa  and  a 
10  bit  exponent  Below  are  shown  the  various  AP  instructions  that  can  be  used  to 
read  and  treat  these  fields  as  three  variable  length  integers  Note  that  some  of 
the  bits  are  inaccessible  in  this  manner. 
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processor,  we  average  somewhat  less  than  50%  This  is  primarily  due  to  the 
extensive  logic  and  indexing  necessary  to  handle  the  general  case  V.uch  of  this 
code  is  overlapped  with  the  actual  floating  point  computation  (since  the  logic 
consists  of  integer  operations)  With  a  conventional  computer,  the  integer 
operations  would  have  to  be  done  in  serial  with  the  floating  point  operations 
This  group  has  also  implemented  molecular  dynamics  calculations  for  water*2 
and  idonte  Carlo  calculations  for  liquid  ammonia*2  and  for  liquid  methane  **  In 
addition.  Andersen  and  Swope*0*5  have  implemented  molecular  dynamics  code 
for  water  and  atomic  solutes  in  water,  and  Berne  and  co-workers*®  have  carried 
out  Vonte  Carlo  calculations  both  on  an  AP-1203  Dammkoehler  and  co-workers*7 
are  implementing  molecular  mechanics  code  on  multiple  cs?;*6  array  proces¬ 
sors  Also  of  related  interest  are  the  crystallographic  work  of  Furey  end  co¬ 
workers*6  and  the  simulations  of  plasma  dynamics  by  the  UCLA  group,40  both  on 
FP5  AiM209's 

Newton  contains  seven  separate  host  programs  which  operate  cooperatively 
as  shown  in  Fig  5  The  first  of  these  is  APCOU  which  is  an  interactive  command 
interpreter  It  allows  the  user  to  type  in  various  commands,  validates  the  com¬ 
mands  and.  where  appropriate,  instructs  the  next  link  m  the  chain  a?. TJX,  to 
perform  the  operation  requested  A?R’JN  is  the  program  responsible  for  running 
the  AP,  and  its  only  function  is  starting  and  stopping  the  AP  and  emptying  the 
data  collection  buffers  while  the  AP  is  running  It  uses  two  other  programs 
APLOA D  to  load  the  AP  initially  with  constants  and  parameters,  and  AP7A3-2S  to 
calculate  the  force  lookup  tables  for  the  intermolecular  interactions  whenever  a 
new  system  is  simulated  DRAW  is  an  optional  process  that  can  be  used  to  display 
the  dynamics  on  the  Picture  System  as  they  are  being  calculated  on  the  AP 
MOVIE  can  be  used  to  display  the  previously  calculated  dynamics  of  a  system 
using  coordinates  saved  m  a  file  by  APRUN  It  also  allows  color  movies  to  be  made 
from  these  files  DODATA  is  a  program  which  allows  the  computation  in  back¬ 
ground  mode  of  a  series  of  runs  to  collect  data  over  an  ensemble  of  initial  start¬ 
ing  configurations  and  momenta  It  is  interruptible  to  allow  users  to  carry  out 
program  testing  and  development  and  systems  tasks  such  as  disk  backup 
without  interfering  with  the  runs 

AH  of  the  actual  code  involved  in  molecular  dynamics  calculation  has  been 
written  in  AP-120B  assembly  language  and  split  into  independent  modules,  as 
shown  in  Fig  6  Other  AP  modules  not  described  here  provide  for  other  types  of 
molecular  mechanics  and  for  the  calculation  of  phenomena  derived  Trom  molec¬ 
ular  dynamics,  such  as  infrared41'42  Raman.  4^'*°  and  electronic46'47  spectra  and 
thermodynamic  quantities  46  The  modules  are  then  linked  together  using  a  sim¬ 
ple  vector  function  chainer  program  that  loops  over  the  routines  to  perform  the 
number  of  integration  steps  requested  Data  is  buffered  inside  the  AP-1203  and 
DMA  ed  out  asynchronously  when  the  buffers  are  filled,  while  the  AP  is  running 
The  buffering  mechanism  will  be  examined  later  in  more  detail 

The  remainder  of  this  section  concentrates  on  our  approach  to  molecular 
mechanics  The  implementation  of  molecular  dynamics  consists  of  two  basic 
■teps.  force  evaluation  followed  by  numerical  integration,  ft  is  these  two  func¬ 
tions  which  are  performed  by  the  AP  modules  shown  in  Fig  6  arid  the  methodol¬ 
ogy  behind  their  operation  and  design  which  will  be  discussed  We  examine  the 
types  of  lists  and  indices  needed  for  a  general  purpose  molecular  mechanics 
package  We  also  examine  various  types  of  boundary  conditions  which  are 
Important  in  many  applications,  especially  those  with  long  range  forces 


computer  TTielr  Intercommunication  and  Interface  to  the  AP  are  shown  On  the 
right  Is  an  optional  Interface  from  the  A P  to  another  host  computer,  a  DKC  PDF 
n '34,  to  which  an  Evans  and  Sutherland  Picture  System  Is  attached  The 
modules  on  the  right  are  C  program  packages  that  run  on  the  PDP  11/14 
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How  to  Describe  an  Atom 

In  addition  to  initial  coordinates  and  momenta,  other  information  is  needed 
to  create  the  wide  variety  of  lists  necessary  to  look  up  masses,  positions,  and 
forces  among  atoms  These  can  all  be  derived  from  a  basic  set  of  information. 

■tract  [ 


char 

fla£a, 

char 

type. 

int 

pa-enC, 

tnt 

param. 

1  Peru, 

where  the  above  is  the  C  language49  template  for  a  structure1  describing  the 
atom  parts  The  first  item  in  the  structure  is  the  atom  flags  These  flags  car.  be 
used  to  enable  or  disable  particular  features  For  example,  one  flag  bit  can  be 
■et  to  fix  the  atom  in  space  so  that  it  cannot  move  (although  it  will  still  exert 
forces  on  neighboring  atoms)  Another  is  used  to  specify  that  the  atom  is  the 
start  of  a  molecule  If  this  flag  is  set  then  the  param  pari  contains  the  total 
number  of  atoms  that  follow  which  are  in  the  same  molecule  Another  flag  b.t 
can  be  set  to  indicate  that  the  atom  is  part  of  a  ring  structure  In  this  case  the 
param  element  contains  the  other  parent  of  this  atom  necessary  to  complete 
the  ring  structure  Another  flag  is  used  to  specify  if  the  atom  is  to  be  drawn  in 
the  picture  display  or  not  Currently  these  four  flags  are  the  only  ones  used, 
although  B  such  flags  are  possible  for  future  uses 

Another  part  of  the  structure  is  the  type  of  the  atom,  a  number  between  C 
and  255  For  example  for  water  the  atom  types  are  hydrogen  and  oxygen  In 
organic  molecules  or  ionic  solutions  it  is  often  necessary  to  distinguish  between 
different  types  of  multivalent  species  (such  as  carbon)  or  different  ionic  state? 
and  thus  each  different  chemical  state  of  a  particular  elements  will  have  a 
different  type  number 

The  parent  element  specifies  to  which  atom  the  present  atom  is  connected 
Using  this  tree  Unking  mechanism  the  entire  boDd  structure  of  any  non-ring 
molecule  can  be  determined  With  the  addition  of  the  ring  flag  and  the  extra 
link  provided  by  the  param  word,  any  reasonable  chemical  structure  can  be 
handled 

The  ato mflag.  type  and  parent  are  stored  as  the  exponent,  high  mantissa 
and  low  mantissa,  respectively,  of  an  AP-ia>B  writable  table  memory  location 
The  atom  param  word  occupies  the  low  mantissa  of  a  main  data  memory  word 
whose  other  fields  are  currently  unused 

Solely  from  this  simple  information  all  the  other  lists  needed  by  the  AD-i203 
modules  described  below  can  be  generated  The  atom  parts,  atomic  coordi¬ 
nates.  and  momenta  are  the  only  pieces  of  information  kept  in  Newton  fill  fijes 
which  are  used  to  start  or  reload  a  Newton  run 

Boundary  Conditions 

Two  of  the  AP  modules  (the  intermolecular  force  evaluater  and  the  integra¬ 
tor)  depend  on  the  type  of  boundary  conditions  being  used  Currently  we  have 
programs  that  allow  the  use  of  four  types  of  boundary  conditions  soft  walls, 
hard  walls,  minimum  image  cubic  periodic,  and  minimum  image  truncated 
octahedral  periodic  00  Other  possible  boundary  schemes7  include  spherical  and 
periodic  boundary  conditions  using  Ewald  sums. 

Cubic  soft  walls  are  the  simplest  No  Imaging  is  done,  and  thus  particles 
feel  only  the  forces  of  the  other  particles  In  the  cube  When  a  particle 
approaches  a  wall,  a  soft  spring  force  pushes  the  particle  back  into  the  cube 
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Ttus  type  of  boundary  condition  is  useful  for  studying  small  clusters  or  droplets 
of  particles  The  disadvantage  is  that  a  high  fraction  of  the  particles  can  be  on 
the  surface  of  the  droplet  or  cluster  and,  in  many  situations,  these  surface 
effects  can  be  important  Consequently,  a  larger  number  of  particles  are 
needed  to  study  bulk  phenomena  In  addition  collisions  with  a  wall  can  give  a 
molecule  a  large,  artificial,  angular  momentum  as  it  is  shoved  back  toward  the 
cluster 

Hard  walls  are  specularly  reflecting  and  cause  problems  with  integration  as 
the  velocity  of  the  particle  normal  to  the  wall  reverses  itself  in  one  integration 
time  step  Such  a  discontinuous  change  can  cause  integration  algorithms  to 
"blow  up"  This  can  be  avoided  by  altering  the  algorithm  so  as  to  alter  the  past 
time  history  of  the  particle  (as  far  back  as  necessitated  by  tbe  algorithm)  when 
it  strikes  the  wall  to  become  that  of  a  parUcle  having  entered  the  box  with  the 
reversed  normal  velocity.  The  problems  with  surface  effects  still  remain 

Periodic  boundary  conditions7  are  commonly  used  to  reduce  surface  effects 
for  the  simulation  of  bulk  matter  The  simplest  is  a  cubic  minimum  image  In 
this  scheme  the  system  of  particles  resides  in  a  central  cube  which  is  sur¬ 
rounded  by  exact  replicas  of  this  central  cube  on  all  sides,  edges,  and  corners 
Particles  interact  only  with  the  closest  image  of  any  other  particle  In  all  cubic 
periodic  boundary  algorithms,  when  a  particle  leaves  the  central  box  it  is 
replaced  by  one  of  its  images  entering  from  the  opposite  side 

A  truncated  octahedral  boundary  condition50  is  similar  except  that  the  unit 
cell  is  a  truncated  octahedron  which  more  closely  resembles  a  sphere  Thus  is 
important  for  minimum  image  boundary  conditions  as  the  forces  must  be 
smoothly  feathered  to  zero  at  the  radius  of  the  inscribed  sphere  of  the  unit  cell 
to  avoid  abrupt  changes  in  force  and  loss  of  energy  conservation  For  a  cube, 
48%  of  the  volume  of  the  cube  lies  outside  the  inscribed  sphere  in  the  corners  of 
the  cube  and  particles  in  this  volume  don't  contribute  to  the  forces  on  the  test 
atom  For  the  truncated  octahedral  geometry,  only  4  5%  of  the  unit  cell  volume 
lies  outside  the  inscribed  sphere,  and  thus  more  of  the  dynamics  calculation  is 
effectively  used  In  addition  the  excluded  volume  is  more  evenly  distributed  in 
angle  than  for  the  cube,  and  the  isotropy  of  space  is  thus  less  distorted  There 
exists  an  easy  way90  to  code  algorithms  to  implement  the  truncated  octahedron 
The  number  of  possible  space  filling  solid  tessellations  is  small  Out  of  the  regu¬ 
lar  and  Archimedean  polyhedra  there  are  only  five  which  are  space  filling  the 
cube,  triangular  prism,  hexagonal  pnsm,  rhombic  dodecahedron,  and  truncated 
octahedron,2' 51  and  thus  the  natural  alternative  to  the  truncated  octahedron 
would  be  the  rhombic  dodecahedron 

Inlermolecular  Force  Evaluation 

The  first  module  of  AP  code  is  the  intermolecular  force  evaluater,  used  to 
compute  nonbonded  forces,  i  e  those  between  atoms  on  different  molecules,  or 
separated  by  so  many  bonds  in  a  single  molecule  as  to  be  considered  indepen¬ 
dent  As  pointed  out  above,  It  comes  at  present  in  four  flavors  soft  and  hard 
walls,  cubic  periodic,  and  truncated  octahedral  periodic  boundary  conditions 
We  approximate  all  intermolecular  forces  as  purely  pairwise  additive,  thus 

MDknljd 

^T^).  a) 

w 

to  which  r^ig  the  distance  between  the  t  th  and  j  th  atoms  The  potential  func¬ 
tion  V{t^)  depends  solely  on  the  chemical  type  of  the  atoms  involved 
Currently,  intermolecular  forces  are  evaluated  by  looping  over  all  tbe  possible 
pairs  of  atoms  in  the  system  simultaneously  calculating  the  force  for  both 
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members  of  tbe  pair.  This  requires  logic  that  allows  the  intermolecular  force 
evaluater  to  skip  over  all  the  pairs  of  atoms  whose  forces  are  to  be  calculated  by 
one  of  the  other  bonded  force  routines 

To  skip  over  the  appropriate  bonded  interactions  the  parent  (and.  if  the 
ring  flag  is  set,  the  param  word)  is  used  to  determine  bonding  up  to  and  includ¬ 
ing  four  body  interactions  which  causes  the  intermoleculer  force  evaluater  to 
■kip  over  these  interactions  While  this  involves  extensive  integer  arithmetic 
and  logic,  the  overhead  is  inconsequential  as  it  is  overlapped  completely  with 
the  code  to  do  the  minimum  imaging  for  the  periodic  boundary  conditions  Ear¬ 
lier  experimental  versions  of  the  inte molecular  module  did  not  carry  out  this 
logic  and  instead  calculated  intermolecular  forces  for  all  pairs  of  atoms  whether - 
bonded  or  not  The  bonded  force  evaluaters  then  simply  subtracted  out  these 
erroneous  intermolecular  forces  when  they  calculated  the  intramolecular 
forces.  This  addition  and  subtraction  caused  disastrous  results  due  to  numeri¬ 
cal  round  off  caused  by  adding  the  relatively  small  intramolecular  forces  to  the 
large  erroneous  intermolecular  forces  that  had  not  yet  been  subtracted  out 

Andersen30  has  pointed  out  that  when  using  a  smoothing  function  for  poten¬ 
tial  energy,  the  correct  force  evaluation  involves  calculating  both  the  force  and 
potential  energy,  as  can  be  seen  from. 

K(r)  =  V(r)S(r)  (2) 

F  =  -  =  -  -^-S(r)  +  V{r) 

ar  dr 

where  V^(r)  is  the  V(r)  potential  smoothed  Lo  zero  by  the  smoothing  function 
5(r). 

Currently  all  intermolecular  forces  are  calculated  by  table  look  up  This 
involves  allocating  most  of  main  data  memory  to  the  force  look  up  tables  Linear 
interpolation  cf  these  grid  points  is  used  to  give  the  actual  forces  At  least  in 
principle,  as  has  been  pointed  out  by  Andersen  and  co-workers,52  this  scheme 
has  a  flaw  when  used  with  the  Verlet  integration  algorithm  due  to  the  infinite 
second  derivative  of  the  force  at  the  boundary  between  the  linear  segments 
Andersen  uses  a  better  scheme  employing  a  polynomial  to  fit  fixed  length  seg¬ 
ments  of  the  potential  curve  with  each  polynomial  joining  smoothly  and  continu¬ 
ously  out  to  several  derivatives  at  the  end  points  of  the  adjoining  segments  This 
involves  less  data  storage  to  calculate  the  intermolecular  forces  as  only  the 
polynomial  coefficients  need  be  stored  In  addition,  since  polynomials  are  used, 
the  potential  energy  as  well  as  the  forces  can  be  calculated  with  little  extra 
effort  from  the  same  set  of  coefficients 

For  large  systems,  most  of  the  computational  time  is  spent  in  intermolecu¬ 
lar  force  evaluation.  6ince  when  done  on  a  pairwise  basis  for  the  entire  system  of 
ff  atoms  it  becomes  a  calculation  proportional  to  Af8,  while  the  bonded, 
intramolecular  calculations  only  scale  as  Af  A  possible  alternative  is  to  use 
neighbor  lists3-53  which  are  only  updated  after  several  time  steps  so  that  only 
the  atoms  which  are  close  to  the  atom  in  question  are  scanned  lo  calculate  the 
Intermolecular  forces  However,  neighbor  lists  create  a  considerable  storage 
problem  In  that  each  atom  must  have  a  list  of  all  the  other  atoms  which  are  near 
It.  For  a  relatively  large  system  this  list  could  easily  exceed  the  amount  of  main 
data  memory  available  For  small  systems,  the  effort  involved  m  updating  and 
indexing  the  interactions  from  such  a  list  may  exceed  the  effort  to  do  all  the 
pair-wise  calculations. 
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Two  Body  Module 

The  two  body  module  calculates  all  simple  bonded  forces  w  diatomics  such 
as  CO  and  Nt  The  force  between  two  bonded  atoms  t  and  j  with  coordinates  r, 
and  ij .  respectively,  is 

Ft  =  -F,  =  V  =  ?  V  =  -  .  (4) 

“rv> 

In  which  Fx  and  Fj  are  the  vector  forces  on  the  » th  and  j  th  atoms,  respectively, 
which  are  separated  by  the  distance  rv  =  ir,-r^|  is  the  gradient  with 
respect  to  the  Cartesian  coordinates  of  atom  i  as  expressed  in  Eq  (Al)  of  the 
Appendix 

The  two  body  program  uses  a  list  as  shown  in  Fig .7.  Hie  low  mantissa  and 
exponent  of  the  first  word  in  the  list  are  used  to  index  the  two  atoms  involved 
Hie  atom  number  of  the  first  atom  (stored  in  the  low  mantissa)  is  subtracted 
from  the  atom  number  of  the  second  before  it  is  stored  in  the  exponent  field 
This  allows  us  to  use  the  fields,  such  as  the  exponent,  to  index  an  atom  number 
which  in  principle  can  be  much  larger  in  magnitude  than  the  bit  field  could  nor¬ 
mally  handle  Since  atoms  that  are  bonded  to  each  other  tend  also  to  be  close 
to  one  another,  the  atom  numbers  are  very  close  in  magnitude,  and  thus  the 
smaller  bit  fields  are  big  enough  to  allow  this  relative  indexing  of  atoms  Code 
common  to  all  two  body  force  evaluations  is  used  to  fetch  the  atomic  coordi¬ 
nates  and  calculate  the  internuclear  vector.  The  potential  index,  or  switch 
parameter,  is  then  used  to  pick  among  a  variety  of  force  calculators  (such  as 
harmonic.  Morse.  etc.)  that  will  calculate,  given  the  internuclear  vector,  the 
scalar  force  along  the  bond  The  low  mantissa  of  the  second  word  in  the  list  is 
the  address  of  the  force  constants  for  the  force  evaluator  to  use  In  this  way.  for 
example,  one  routine  can  be  used  to  evaluate  all  harmonic  forces.  The  common 
main  line  code  then  decomposes  the  force  along  the  space  fixed  x.  v-  and  z 
axes  This  module  is  only  used  for  diatomics  For  more  complex  molecules  it  is 
more  efficient  to  have  the  three  body  module  also  calculate  the  two  body  forces 


three  Body  Module 

Three  body  interactions  are  those  whose  calculation  depends  on  the  coordi¬ 
nates  of  three  particles,  r,.  r j .  and  rt  An  example  of  this  is  a  force  due  to  bond 
angle  bending  which  requires  the  three  atoms  involved  to  be  specified  in  order 
to  calculate  the  bond  angle,  d.  The  three  body  module  uses  a  list  similar  to  the 
two  body  module,  as  Bhown  in  Fig  7.  The  low  mantissa  of  the  first  word  contains 
the  index  to  the  middle  atom  in  the  three  body  interaction  The  exponents  of 
the  first  and  second  words  contain  the  atom  Dumbers  of  the  other  two  atoms 
after  subtracting  the  atom  number  of  the  middle  atom  The  high  mantissa  of  the 
first  word  selects  which  three  body  force  evaluation  routine  is  to  be  used 
Currently  we  have  two  such  evaluaters.  a  complex  one54  for  water  molecules  and 
a  simpler  harmonic  one, 

K(dr.,<5r»,dd)  =  Aq(5r.)B  +  +  A* (W)2  ,  (5) 

which  is  written  In  terms  of  the  bond  vectors  r,  and  r*  where 

tfr.  =  |rt  ~r}  j  -r*  =  |r.  |  -r»  (6) 

6**  -  I **  -*)  I  ~ri  =  I**  |  ~rg  ,  (7) 

and  id  =  d— d* ,  In  which  rj,  rf  and  d*  are  the  equilibrium  bond  distances  and 
angles,  respectively,  of  the  potential  V.  The  two  body  part  of  the  potential  is 
shown  in  the  above  expression  as  it  is  also  calculated  in  this  module  as  explained 
above  The  water5*  force  evaluator  has  various  higher  order  terms  among  6rm . 
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Figure  7.  Format  of  interaction  lists  The  top  pane]  shows  the  format  of  the  two 
body  interaction  list  as  stored  in  A P  main  data  memory  The  low  mantissa  of  the 
first  word  holds  the  i  th  atom  number  times  three  (for  faster  indexing  into  the 
three  dimensional  arrays)  The  high  mantissa  is  an  B  bit  integer  specifying 
which  force  calculator  to  use  (harmonic,  Korse,  etc).  The  exponent  field  con¬ 
tains  the  difference  between  the  j  th  and  <  th  atom  numbers  This  is  done  to 
allow  more  dynamic  range  in  the  10  bit  field  The  second  word  contains  the 
address  of  the  force  constant  The  switch  parameter  or  potential  index  is  used 
to  specify  which  force  evaluation  routine  is  to  be  used  for  this  interaction  The 
three  body  list,  in  the  middle  panel,  is  similar  except  that  the  middle  atom  is 
used  to  Index  the  other  two  The  four  body  list,  as  shown  in  the  bottom  panel,  is 
also  similar  except  for  an  constant  offset  which  is  subtraced  from  the  force  con¬ 
stant  pointer  to  allow  more  dynamic  range. 
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6rt ,  and  66  in  addition  to  the  terms  shown  m  Eq  (5)  The  low  mantissa  of  the 
second  word  contains  the  address  of  the  force  constants  (*o.  *i-  and  *2)  used  by 
the  force  evaluators 

Common  code  is  provided  by  the  three  body  module  to  calculate  the  two 
internuclear  vectors  and  the  bond  angle  before  calling  the  appropriate  three 
body  force  evaluator.  The  evaluator  returns  scalar  forces  along  the  two  inter¬ 
nuclear  vectors  and  a  force  associated  with  the  bond  angle,  and  the  three  body 
module  resolves  these  forces  into  Cartesian  forces  on  the  three  atoms  involved 
The  forces  on  the  atoms  are  given  by 

*  =  -V.  v  .  -  VV.r.  *  !£*.,.]  =  -  j*.V  ♦  ££».«]  (6) 

r.  •-V.v.  -  -  -jf.V  (10) 

where  we  have  used  Eqs  (A9)  and  (A10)  of  the  Appendix  to  evaluate  the  chain 
rule  gradients  for  the  tensor  vector  product  appearing  in  Eqs  (8)-(i0).  Simi¬ 
larly.  using  the  results  of  the  Appendix,  the  gradients  involving  the  bond  angle  i5 
in  Eqs  (8)-(10)  can  be  expanded  in  terms  of  gradients  involving  the  bond  vectors 

r,  and  r* . 

f,6  =  =  V.6  (ll) 

Wjif  =  Vm*Vjra  +  =  -?.6-V*6  (12) 

?*6  =  ¥*6V*r*  (13) 

Since  the  bond  angle  can  be  written  in  terms  of  the  dot  product  of  the  two  bond 
vectors. 


COST?  = 


r.  r> 


the  terms  in  Eqs  (11)-(13)  can  be  evaluated  as 

W  Jt  —  «F 


Using  Eqs  (14)  and  (15). 


▼,  cost!  = 


where 


,cosd  . 


▼.  (rB  r*  )  -  r.  r* Wm  (r.r*  ) 


▼.(*■.  r*)  =  r* 


Wm(rmrt)  =  ~r.  . 

Tm 

Substituting  Eqs  (17)  and  (18)  into  Eq  (16)  yields 


V,cos6  =  - 

fl 

Now.  by  using  Eq  (14)  with  Eq  (19)  we  have 


%r*r*  -  ~(rm  r*)r. 
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W.COSd  =  — - 


record 


(23) 


A  similar  derivation  for  the  gradient  with  respect  to  the  other  bond  vector  rt 
yields 


V^cosiS  = 


r#  cosiS 


r«r» 


Finally,  substituting  Eqs  (20)  and  (21)  into  Eq  (15)  we  have 


= 


-1 

r» 

racosd 

sim5 

Tart> 

r.e 

r. 

r„c  osd 

Bind 

r,rt 

rS 

(21) 


(22) 

(23; 


Thus.  using  Eqs  (8)-(13)  and  Eqs  (22)  and  (23)  the  forces  can  be  appropriately 
resolved  onto  the  three  atoms  involved  These  results  are  the  same  as  those 
arrived  at  using  the  Eliaschench  and  Wilson  e-vector  method66’66  to  evaluate  the 
elements  of  the  B  matrix  used  in  normal  mode  vibrational  analysis  to  relate  the 
internal  and  Cartesian  coordinates  through  a  Taylor's  expansion 


Fbur  Body  Module 

A  four  body  interaction  requires  the  knowledge  of  the  positions  of  four  par- 
ticles  to  calculate  the  force  The  two  most  common  examples  are  torsional 
forces  and  out-of-plane  bending  forces5656  The  four  body  module  for  torsional 
forces  uses  a  list  as  shown  in  Fig  7  The  atom  numbers  of  the  two  inner  atom? 
are  stored  in  the  low  mantissa  field  of  the  first  and  second  words,  with  the  j  th 
atom  number  multiplied  by  three  for  indexing  convenience  The  inner  atom 
index  numbers  are  first  subtracted  from  the  closer  outer  atom  number  and  then 
stored  in  the  exponent  field  of  the  two  words  The  high  mantissa  of  the  first 
word  is  used  to  store  the  potential  index  value  to  select  a  particular  torsional 
force  evaluator  with  the  force  constants  being  indexed  by  the  value  in  the  high 
mantissa  field  of  the  second  word  Since  a  complete  memory  address  cannot  be 
stored  here,  the  value  in  the  high  mantissa  of  the  second  word  is  an  offset  to  the 
base  address  of  the  torsional  force  constants  Currently  there  are  three  types 
of  torsional  force  evaluaters  which  handle  single,  double,  and  triple  bonds 

Common  code  16  used  to  calculate  all  the  internuclear  vectors  and  the  tor¬ 
sional  angle.  A  scalar  force  is  returned  which  is  solely  dependent  on  the  tor¬ 
sional  angle  This  force  is  then  decomposed  onto  the  atoms  as  foCows  If  we 
take  a  four  body  interaction  as  shown  in  top  half  or  Fig  B  each  of  the  four  atoms 
<,  j ,  k .  and  l  has  coordinates  represented  by  the  vectors  rt .  r; .  r* .  and  r(  If  we 
now  define  the  bond  vectors, 

*T  =  1-1-17,  rg  =  r4-r/.  rs  =  r{-r*.  (24) 

and  look  down  the  j-k  bond,  then  we  have  the  representation  shown  in  bottom 
h"lf  of  Fig  8  where  r»  and  rt  are  the  projection  of  rt  and  rs  into  a  plane  and  the 
tf  torsional  angle  is  formed  between  them  We  can  now  write 

r«  =  r,  Xri  r.  =  |r.  |  .  (25) 

r»  =  rjXr*  r*  =  |rv  1  .  (26) 

and 


cosjp  = 


(r,  X  r2)  (rs  X  r2) 


r.  r, 

rmU 


(27) 


|rjXr2|  |rs X rs | 
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The  force  on  each  atom,  using  Eqs  (A9)  and  (AID)  from  the  Appendix,  is  given  by 


-?,v=  -V.l'V.r,  =  -V,V  (26) 

-V,vvjri  -  V2l'V2  =  ^V  +  VzV  (25) 

-¥*  V  =  ~VzVVtr2  -  VjKV*  rs  =  -VZV  +  VSK  (30) 

-f(V=  -VjVV.ra  =  -?3V.  (3:) 


Since  the  potentials  used  are  solely  functions  of  the  torsional  angle. 


V  =  V(cosit)  ,  (30) 

each  of  the  terms  in  Eqs  (28)-(3l)  can  be  evaluated  as 

vv=d^coS£]_  {=1.2.3.  (33) 

*  cfcosp  * 

However,  following  the  format  of  Eqs  (A4)-(A7)  of  the  Appendix,  we  can  expand 
the  gradients  of  cosjf  in  terms  of  the  projections  r»  and  to  give 


V{cos^>  =  V.  cos^r  V<r.  +?.ces^V<r*  .  {=1,2.3  (34) 

From  Eqs  (20)  and  (21)  we  know  the  gradients  of  the  torsional  angle  are  given 
by 


V»cos?  = 


r„cos? 


(3 


V*cosjf  = 


rbCOS? 


(36) 


Tf 

Using  Eqs  (35)  and  (36)  with  the  vector  tensor  products  of  Eq  (34)  which  are 
given  by  Eqs  (Al3)-(A15)  of  the  Appendix,  we  have 

V,V  =  r*X  \~^~ 


r0eosjr 


?EV  = 


r* 

r.cosp 

X  r  4- 

r„ 

r„  cos? 

r«r* 

r! 

A  Ij  4- 

U* 

r3 


*3F  =  rE  X 


r^cos? 


r.r* 


(37) 

(36) 

(39) 


Using  Eqs  (37)-(39)  with  Eqs  (28)-(3l)  the  forces  on  all  four  atoms  can  be 
resolved 


Integration  Module 

The  second  step  in  molecular  dynamics  is  the  numerical  integration  of 
Newton's  Second  Law  for  each  atom  i, 


-  dtr' 

m,  "*  dt/B  ‘ 


(40) 


where  m*  is  the  mass  of  the  x  th  particle,  St  its  acceleration.  r»  its  position,  and 
t  time.  As  pointed  out  above,  the  integration  module  depends  on  the  type  of 
boundary  condition  in  use,  aB  this  module  also  applies  any  position  changes 
necessary  to  keep  the  atoms  In  the  unit  cell  (in  the  case  of  periodic  boundaries) 
or  applies  any  restoring  forces  necessary  (in  the  case  of  soft  walls)  The  integra¬ 
tion  algorithm  we  use  in  each  of  these  modules  is  the  same,  a  version  of  the  Ver- 
let  algorithm  as  discussed  by  Bee  man57  with  further  modifications  by  Ander¬ 
sen.®  In  our  implementation,  the  vector  difference  in  positions  d,(f  +  1)  at  time 


k 


step  t  +  3  is  calculated  from  the  previous  difference  d,(0  tad  the  force. 

d +  =  4(0  * 

to, 

In  which  t  indicated  the  time  step  t .  and  A  is  the  size  of  the  time  step  Next  the 
new  positions  are  calculated  from  the  new  difference  in  positions 

rx[t  +  1)  =  rt(<)  +  d,(f +  l)  («) 

The  improvement  made  by  Andersen  is  that  now  only  the  difference  in  position  is 
stored  rather  than  the  velocity  Therefore  there  is  less  round  off  error  as  the 
numbers  being  added  in  each  stage  of  the  calculation  are  closer  in  magnitude 
A  crude  forward  difference  velocity  can  easily  be  obta.ned  by  dividing  the 
difference  in  position  by  the  time  step,  or,  if  desired,  more  accurate  velocities 
can  be  calculated  67 

Beeman57  shows  that  higher  order  integration  techniques  tend  not  to  be  as 
stable  as  the  simple  Verlet  algorithm  when  larger  time  steps  are  used  He  also 
shows  that  the  Verlet  algorithm  conserves  energy  as  well  as  other  integrators 
tested  in  the  large  time  step  Limit  The  advantage  of  using  this  simple  techmquc- 
with  the  array  processor  is  that  only  a  minimal  amount  of  memory  is  set  aside 
for  the  storage  of  positions  and  the  past  time  history  of  the  system  as  compared 
with  higher  order  integration  methods 

Data  Collection  Modules 

The  data  collection  routine  used  for  a  particular  application  is  often  very 
specialized  and  thus  to  analyze  different  properties  different  AP  modules  must 
be  written  However,  IDO'JT,  the  mechanism  for  buffering  the  data  out  to  be  the 
host  computer,  is  common  to  all  routines  once  the  property  that  is  to  collected 
is  calculated  In  many  cases  this  property  can  be  calculated  in  a  vector  function 
program 

Various  data  collection  modules  exist  for  use  with  Newton  An  example  is 
the  module  that  calculates  and  collects  dipole  moments  and  polarizabilities  for 
the  system  as  a  function  of  time,  in  order  to  compute  by  linear  response  theory 
the  infrared41’42  and  Raman 4^~45  spectra  These  modules  operate  in  a  very  simi¬ 
lar  manner  to  the  force  evaluation  modules,  in  that  for  the  most  part  they  use  a 
hst  to  index  the  atoms  involved  in  the  data  calculation  and  the  relevant  parame¬ 
ters  (such  as  static  and  derivative  dipole  moments  for  molecules)  IDO'JT  saves 
the  data  in  two  Internal  buffers  m  main  data  memory  and  then  DMA's  the  data  out 
under  interrupt  control  using  a  double  buffering  technique  described  in  the  next 
section  In  general,  the  actual  data  is  all  that  is  saved  from  an  individual  Newton 
run  Normally  trajectories  are  not  saved  in  that  it  is  easier  and  faster  to  do  the 
run  over  rather  than  saving  and  storing  the  details  of  the  run 

TV  RESULTS  AND  ANALYSIS 

As  pointed  out  in  section  D.  the  main  advantage  of  an  array  processor  is 
high  computational  speed  at  a  relatively  low  price  In  the  way  of  benchmarks, 
we  have  obtained  the  following  results  presented  in  Table  11  Column  one  shows 
the  various  computers  on  which  the  benchmarks  are  taken  Columns  two,  three, 
and  four  compare  the  speed  of  the  AP-120B  for  three  different  molecular  mechan¬ 
ics  packages  The  one  presented  in  column  two  is  a  direct  C  translation  of  our 
AP  molecular  dynamics  package  run  on  a  VAX  lirroo  and  aVAX  n  nyo  without  float¬ 
ing  point  accelerators  (FPA's)  Column  three  compares  our  AP  version  against  an 
optimized  Fortran  version  due  to  Hagler5®  and  run  on  a  VAX  iimso  with  a  floating 
point  accelerator  Column  four  compares  some  speeds  obtained  for  lfonte  Carlo 
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calculations  from  Chester,  et  a L  59 

Table  H.  Compar-jon  a?  S_rr:  jle: ion  ?..-r.es  on  Various  CoT.pu.rj 


Corr.p  jler 

Mo  ec  Sji' 
fb  C  m  c  oh  VAX' 

■  Dyri£*T_c3 

•  •  r?A  ta  VAX' 

Modi*  Ct-.o 

AP-1203 

i 

1 

1 

VAX  11/780 

BO 

35 

VAX  11/7S0 

160 

PD?  11 

300 

Pricit  400 

65C 

IBM  37G  -'166 

16 

CDC  7600 

0  47 

As  car.  be  seer  from  the  tabie  our  molecular  mechanics  package  is  appro:.- 
mate'.y  thirty-five  times  faster  than  a  VAX  n/780  with  a  floating  point  accelerator 
and  an  optimized  Fortran  compiler  Thus  a  simulation  that  can  be  run  in  a  wee  k 
and  a  half  on  the  A?- 1203  would  take  a  year  on  a  VAX  even  if  the  vax  were  totally 
dedicated  to  that  calculation 

Although  the  AP-1203  has  proven  to  be  a  very  fast  and  economical  machr.-: 
for  molecular  mechanics  and  allows  us  to  simulate  systems  which  otherwise 
would  not  be  feasible  it  is  far  from  idea!  One  concern  is  the  word  length  of  the 
floating  point  numbers  For  example,  for  many  quantum  mechanical  calcula¬ 
tions  32  bit  floating  point  representations  are  inadequate  whereas  64  bit  pre-c  - 
sion  suffices  The  question  arises  as  to  the  adequacy  of  the  38  bits  of  the  AP  and 
If  particular  sections  of  algorithms  can  be  painstakingly  coded  in  the  forced 
double  precision  possible  on  the  AP  to  give  adequate  results  While  most  molec¬ 
ular  mechanics  is  certainly  adequately  handled  with  36  bits,  one  might  do  the 
integration  module  for  molecular  dynamics  m  double  precision  and  leave  the 
more  computationally  intensive  force  evaluation  in  single  precision 

Another  feature  missing  in  the  A°-1203  is  the  ability  to  DMA  data  out  of  main 
data  memory  to  the  host  processor  preserving  the  full  3B  bits  of  precision  This 
impacts  the  ability  to  do  three  things  The  first  is  that  we  would  like  to  be  able 
to  rapidly  store  intermediate  results  on  disk  so  that  the  runs  could  be  stopped 
and  started  at  a  later  time  without  loss  of  precision  in  the  data  A  rapid  method 
of  writing  such  intermediate  files  would  also  facilitate  periodic  file  dumps  for 
restart  capability  if  the  host  system  were  to  crash  The  second  impacted  area  is 
the  ability  to  rapidly  load  and  retrieve  integers  packed  into  floating  point  words 
The  amount  of  time  spent  in  loading  the  AP  through  the  virtual  front  panel  with 
long  packed  lists,  for  example  which  would  be  necessary  if  neighbor  lists  were 
used  could  be  extraordinarily  burdensome  Thirdly,  this  limitation  interferes 
with  making  the  AP  a  rapidly  sharabie  machine  rather  than  an  exclusive  use  dev¬ 
ice  as  it  is  now  To  make  it  truly  sharabie  all  of  the  machine's  memories  and 
internal  registers  would  have  to  be  rapidly  DMA'ed  out  to  disk  (swapped)  and  this 
is  not  possible  with  the  present  architecture. 

In  future  versions  of  array  processors  we  would  also  like  to  see  a  separate 
integer  memory  that  could  be  accessed  faster  .  Lacking  this  an  instruction 
should  be  added  that  would  allow  the  access  of  all  the  bits  of  an  integer  packed 
into  a  main  data  word  Since  molecular  dynamics  is  not  usually  memory  bound 
and  most  code  is  parallel  Tor  each  of  the  three  coordinate  axes,  more  than  one 
floating  point  adder  and  multiplier  might  be  useful  In  this  way  calculations  for 
each  of  the  three  axes  could  be  carried  oul  in  parallel  rather  than  serially  The 
serial  approach,  however,  is  often  convenient,  especially  giveD  the  fact  that  the 
multiplier  is  a  three  stage  operation,  but  since  the  adder  is  only  two  stages  it 
often  disrupts  any  attempts  to  build  the  physical  symmetry  of  the  problem  into 
the  code  If  such  Improvements  were  to  be  incorporated,  speed  enhancements 
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by  utilizing  fast  main  data  memory  could  more  easily  be  realized 

Another  obvious  misfeature  is  the  absence  of  integer  multiplication  which 
makes  the  addressing  of  multidimensional  arrays  difficult  Additional  memory 
address  registers  so  that  more  than  one  item  could  be  fetched  from  different 
banks  of  main  data  memory  m  the  same  instruction  would  replace  the  necessity 
of  using  the  writable  table  memory  for  selective  data  storage,  thus  allowing  all 
data  to  be  accessible  through  the  DUa  channel,  freeing  table  memory  for  seldom 
changed  constants 

There  is  currently  available  from.  Floating  Point  Systems  a  64  bit  array  pro¬ 
cessor  known  as  the  F3S-J64  which  solves  the  problem  of  numerical  accuracy,  at 
least  for  most  problems  of  chemical  interest  It  is.  however,  no  faster  (in  fact 
somewhat  slower)  than  the  AP-1203  and  quite  expensive  It  also  differs  in  that 
program  memory  and  data  memory  occupy  the  same  space,  and  thus  the  For¬ 
tran  compiler  approach  although  no  faster  id  speed,  is  more  feasible  due  to  the 
abundance  of  memory  now  available  for  its  bulky  code  This  offers  some  advan¬ 
tages.  foe  example,  for  quantum,  applications  as  large  previously  developed  pro¬ 
gram  packages  can  be  rur.  in  the  AP  using  its  Fortran  comp.ler  with  only  the 
Inner  loops  being  optimized  as  hand  coded  AP  routines  Studies  of  usual  quan¬ 
tum  packages  indicate  that  only  500  to  1000  lines  of  code  take  up  most  of  the 
execution  time  90 

Another  improvement  that  would  facilitate  the  use  of  array  processors  for 
computational  chemists,  is  a  good  higher  level  language  compiler  We  would 
suggest  the  language49  C  For  many  reasons  C  is  more  appropriate  than  Fortran 
for  compiling  into  efficient  A?-i203  assembly  language  For  example  C  allows  the 
use  of  pointers  to  arrays  as  an  alternative  to  subscripts  Incrementing  and 
decrementing  a  pointer  and  stepping  along  memory  can  be  much  more 
efficiently  handled  in  Ap-i203  assembly  language  than  adding  a  subscript  offset 
(which  may  need  to  be  decremented  if  it  does  not  start  at  0.  such  as  is  the  case 
with  Fortran)  to  the  array  base  address  In  addition,  C  allows  the  declaration  of 
variables  as  registers,  thus  allowing  the  programmer  to  warn  the  compiler  that  a 
particular  piece  of  data  needs  to  be  kept  in  date  pad  registers  rather  than  writ¬ 
ten  and  reread  from  memory  Ve  are  optimistic  that  a  good  C  compiler  can  be 
written  that  will  produce  AP-1209  assembly  language  code  good  enough  to 
obsolete  the  desire  to  program  in  assembly  language  except  for  extremely  criti¬ 
cal  loops  which  are  executed  too  many  times  to  tolerate  any  inefficiencies  The 
only  such  loop  or  this  type  in  our  molecular  dynamics  code  is  the  intermolecular 
force  evaluation,  which  is  less  than  ten  percent  of  the  total  AP  code 

The  actual  generation  or  such  a  compiler  is  a  difficult  task  Since  the  AP  is 
not  of  the  von  Neumann  architecture,25  there  is  little  expertise  m  this  area  of 
software  design  Ken  Wilson®1  has  suggested  a  Vonte  Carlo  method  of  code  gen¬ 
eration.  whereby  given  certain  rules  and  constraints  the  AP  itself  would  try  to 
optimize  its  gene-aled  assembly  code  using  a  Vonte  Carlo  technique,  varying 
the  code  while  preserving  its  logical  outcome  In  this  case  code,  not  particles, 
would  be  randomly  moved  with  the  overall  length  of  the  code  being  minimized 
His  attempt  at  imp’ementmg  such  an  optimizer  also  starts  with  C  language 
•ource  code 

All  of  the  support  code  for  Newton  run  on  the  host  computer  is  written  in  C 
This  has  allowed  us  to  maintain  a  single  set  of  source  code  files  which  are  shared 
by  many  programmers  and  used  to  simulate  systems  dramatically  different  in 
nature  C  is  much  more  structured  than  Fortran  and  self  documenting  in  many 
cases,  as  it  has  superior  readability  It  discourages  the  use  of  "go  to"  state¬ 
ments  which  have  been  described  as  a  marvelous  way  to  write  impcssib!e-to- 
understand  programs  Although  it  is  a  higher  level  language  it  allows  the 
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programmer  the  freedom  and  degrees  of  manipulation  of  data  found  tn  most 
typical  assembly  languages  R  is  portable  since  there  is  no  built-in  1/0  and  sys¬ 
tem  dependencies  are  only  in  word  lengths  Although  C  is  very  closely  tied  tc 
the  UNIX  operating  system,  there  are  C  compilers  running  on  VAX's  under  VMS  on 
IBli  360's,  on  and  various  other  machines  There  is  a  portable  C  compiler  that  Cc.n 
be  bought  up  on  most  machines  with  just  a  few  months  of  work  The  C  compiler 
Itself  is  written  in  C 

The  advantages  of  the  UNIX  operating  system®  63  should  not  be  cverl ricked 
The  reason  UNIX  and  C  are  so  related  is  that  all  of  the  UNIX  utilities  and  o\er  9.1 
percent  of  the  actual  UNIX  kernel  are  written  in  C  Thus  the  UNIX  operating,  sys¬ 
tem  itself  is  portable  It  is  becoming  a  standard  operating  system  for  a  w.dc 
variety  of  computers,  so  that  we  can  (and  have)  moved  both  the  UNIX  operating 
system  and  Newton  from  one  type  of  processor  to  another  Its  debugging  edit  - 
mg,  and  friendliness  to  the  user  are  superior,  enhancing  programmer  produc¬ 
tivity  and  the  ease  of  making  and  debugging  changes  to  Newton 

The  development  of  Newton  would  have  been  very  difficult  without  the  L  V  \ 
operating  system  environment  The  operating  system  kernel  ar.d  device  drier.; 
for  UNIX  are  written  m  C,  and  thus  are  easily  changed  We  have  hand  tailored  the 
AP-1203  driver  to  meet  our  needs  It  allows  the  host  computer  and  the  A?-i203  to 
operate  asynchronously,  coordinating  efforts  via  interrupts  As  the  A3-1203  fills 
up  its  data  buffers  it  sends  an  interrupt  to  the  host  which  causes  the  host  tc 
empty  the  buffer  from  AP-1203  memory  to  disk  while  the  AP-1203  continues  the  cal¬ 
culation,  filling  up  a  second  buffer  This  means  that  there  is  no  lost  computation 
time  by  the  AP-1203  waiting  for  the  host  to  empty  the  buffer  and  restart  the  cal¬ 
culation  Furthermore  this  procedure  allows  the  host  program  to  be  inactivated 
(and  even  swapped)  without  having  to  loop  just  to  check  on  the  state  of  the  AP 
This  change  has  increased  our  data  throughput  by  over  a  factor  of  ten  and  allows 
the  use  of  the  AP-1203  on  a  timesharing  system  with  minimal  impact  to  other 
users 

UNIX.  a  third  generation  operating  system,  is  easy  to  learn  Our  research 
group  consists  entirely  of  chemists,  most  of  whom  have  little  previous  training  in 
computer  science  and  most  have  had  little  difficulty  in  picking  up  the  necessary 
■kills  to  use  the  operating  system  on  a  sophisticated  level  Since  most  have  lit¬ 
tle  or  no  experience  in  traditional  computer  languages  such  as  Fortran,  it  is 
interesting  to  note  that  they  can  begin  writing  complicated  C  programs  in  much 
less  time  that  it  would  have  taken  to  gain  the  equivalent  abilities  in  Fortran 

V.  CONCLUSION 

As  computational  chemists  search  for  more  computer  power,  others  will 
surely  turn  to  array  processors  as  we  have,  as  they  provide  at  the  moment  by 
far  the  most  computational  power  per  hardware  dollar,  particularly  since  the 
cost  is  low  enough  that  they  can  be  dedicated  full-time  to  a  particular  task  or 
class  of  tasks  While  running  on  a  supercomputer  such  as  a  Cray-!  will  result  m 
more  computation  per  hour  of  processor  use,  it  is  unlikely  to  result  in  as  much 
computation  per  year  The  reason  is  that  the  equivalent  to  24  hours  per  day  of 
dedicated  A?-iao3  time  is,  for  example,^®4'®'®  two  to  eight  hours  per  day  of 
Cray-1  Ume,  a  usage  rate  which  few,  if  any,  research  groups  are  able  to  afford 
over  the  long  run  Even  if  a  group's  budget  were  large  enough  to  annually  pur¬ 
chase  this  much  supercomputer  time,  for  the  same  cost  several  array  proces¬ 
sors  could  be  purchased  each  year. 

While  array  processor  use  is  very  appealing  and  the  reward  can  be  high,  we 
believe  our  effort  in  bringing  up  a  genera!  purpose  program  package  for  molecu¬ 
lar  mechanics  has  also  uncovered  many  of  the  pitfalls  That  we  can  run  in  ten 
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days  problems  which  would  require  a  year  of  dedicated  VAX  urr&j  time  allows  us 
to  handle  problems  in  solution  reaction  and  biomolecular  dynamics  which  would 
not  otherwise  be  feasible  However,  the  price  we  have  paid  is  substantial  tVru!' 
molecular  mechanics  is  straightforward  in  nature,  it  has  taken  over  six  mar- 
years  to  develop  efficient  AP  code  to  carry  out  the  task 

An  important  feature  of  our  code  is  its  modularity  Since  reprogramming  is 
expensive,  we  have  attempted  to  isolate  the  individual  aspects  of  the  calculation 
into  individual  AP  modules  The  generality  of  the  program  package  allows  us  to 
simulate  a  wide  variety  of  systems  using  essentially  the  same  code  Past  work 
includes  the  calculation  from  molecular  dynamics  and  linear  response  theory  of 
infrared,41'*2  Raman*3"**  electronic*6*7  spectra  in  the  gas  phase  and  in  liqu.d 
solution  In  addition,  we  have  computed  the  dynamics  and  rotational  and  vibra¬ 
tional  spectra  of  alkanes  (such  as  methane,  ethane,  cyclohexane  and  their  solu¬ 
tions),  water  (in  both  the  gas  and  liquid  phase  as  well  as  various  N-roers  of  water 
molecules),  and  ions  and  microcrystals  dissolved  in  water,  We  have  computed 
the  transient  Raman  and  electronic  absorption  spectra  during  the  course  of  a 
chemical  reaction  by  computing  the  dynamics  for  the  photodissoeiation  of 
iodine  in  a  solution  of  liquid  xenon **-*7  Other  applications  involve  the  computa¬ 
tion  from  molecular  dynamics  of  thermodynamic  quantities  and  their  quantum 
correction  through  spectral  analysis  of  atomic  velocity  time  histories  43  Newton 
also  incorporates  a  general  set  of  protein  potentials  for  biomolecules  and  is 
currently  being  appL°d  to  the  molecular  mechanics  of  polypeptides  and  mem¬ 
branes  in  collaboration  with  A  Hagler. 
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APPENDIX  CHAIN  RULE  FOR  GRADIENTS 


In  the  calculation  and  decomposition  of  forces  it  is  often  convenient  to 
switch  Trom  Cartesian  to  internal  coordinate  systems  This  presents  difficulties 
as  the  gradients  must  also  undergo  this  transformation  In  this  appendix  we 
present  formula  for  converting  gradients  in  one  frame  of  reference  to  another 

The  gradient  with  respect  to  the  Cartesian  coordinates  of  particle  i , 
r,  =  *»T+ytJ+  *t£  is  the  vector  operator  given  by 


% 


B_ 

A*. 


(A!) 


The  force  on  the  i  th  particle  F,  can  be  expressed  using  the  operator  in  Eq  (Al) 

as 


Ft  =  ~VXV.  <A2) 

where  V  is  the  potential  energy  Commonly,  however,  the  potential  V  is  more 
easily  expressed  as  a  function  of  some  internal  coordinate  r, ,  where  the  internal 
coordinate  is  a  function  of  the  Cartesian  coordinates,  r,  =  r,  (r,)  We  would 
therefore  like  to  convert  the  gradient  in  Eq  (A2)  into  a  gradient  with  respect  to 
the  internal  coordinate  r,  Using  the  chain  rule,  the  following  terms  result 
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+  [  *v  °T»  ,  +  , 
(a*,  8*,  ay.  9z,  a*,  a*. 

Equation  (A3)  can  be  written  in  a  more  compact  form  as 

%K**)  =  V.  t’(r.)  V,r. 

where  V.  V(r.)  is  a  vector  with  components 


=  [t 

and  V,r0  is  a  tensor64  with  components 
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Thus  by  Eqs  (A5)  and  (A6),  the  expression  in  Eq  (A4)  is  actually  the  vector 
matrix  product65 


9x. 

9x„ 

9x. 

0*i 

9y. 

82, 

9  V 

av 

BV 

»y. 

«y» 

8y. 

9x. 

9y„ 

9z. 

ax. 

9y» 

82, 

9z. 

82. 

dz„ 

ax* 

9y» 

9z, 

which  when  expanded  gives  Eq  (A3) 

We  will  now  apply  this  technique  to  the  examples  presented  in  the  text  In 
the  case  of  the  three  body  module,  where  the  internal  coordinate  vectors  r.  and 
r*  as  in  Eqs  (8)  and  (7)  of  the  text,  respectively,  are  given  by  the  difference  in 
position  of  the  i  th  particle  with  respect  to  the  j  tb  particle,  viz  , 

r,  =  r,  -  ri  (AB) 

It  can  be  easily  seen  by  evaluating  the  partial  derivatives  in  Eq  (A6)  using  Eq 
(AB)  that  the  relevant  tensors  are  given  by 

V,rB  =  I  (A9) 

V^r,  =  -1.  (AiO) 

where  I  is  the  unit  tensor 

However,  in  the  four  body  force  decomposition  the  internal  coordinates  r. 
and  r*  are  vector  cross  products  of  the  bond  vectors  r,,  rfc,  and  ry  viz  , 


r.  =  r,Xr£  =  |yj«£-*iy£.  *,x£-x,*£.  x,y£-y)x£J  (All) 

For  this  case  we  must  evaluate  the  partial  derivatives  in  Eq  (A6)  using  Eq  (All) 
except  with  respect  to  the  Cartesian  vector  ^  instead  of  r,  This  gives  for  the 
gradient  expressed  in  Eq  (34)  of  the  text  the  following  result 
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When  Eq  (A12)  is  expanded  it  becomes  obvious  that  it  is  equivalent  to 

Vicosjf  =  rt  X  V.cosp  .  (A13) 

Similarly,  the  following  results  can  be  obtained  for  the  other  gradients 

VBcos^  =  V.cosy  X  r,  +  W^cosjp  X  r3  (A  1  < ) 

Wacosjc  =  X  V* cosf  (A'.5) 


References 

I.  K  R  Wilson,  in  Computer  Networking  and  Chemistry,  ed  P  Lykos  Ameri¬ 
can  Chemicad  Society,  Washington  1975  p  17 

2  K  R  Wilson,  in  Minicomputers  and  Large  ScaLe  Compulations,  ed  P  Lykos 
American  Chemical  Society,  Washington,  1977,  p  147. 

3  D  Fincham  and  B  J  Ralston,  Cbmput  Ffiys  Comm  .  23,  127  (1981) 

4  D  Pmcham,  Compul  Phys  Comm  .  21,  247  (1BB0). 

5  B  J  Alder.  Arm  Rev  Phys  Chem  ,  24  325  (1973) 

8  W.  W  Wood  and  J  J  Erpenbeck,  Am  Rev  Phys  Chem  ,  27.  319  (1976) 

7.  J  P.  Valleau  and  S  G  Whittington,  in  Sfaftshca!  Mechanics  Part  A  Equil- 

brium  Techniques,  ed  B  J  Berne.  Plenum,  New  York.  1977 

8  1  R  McDonald,  in  Microscopic  Structure  and  Dynamics  of  Liquids,  ed  J 
Dupuy  and  A  J  Dianoux.  Plenum.  New  York.  1977. 

9  K  Binder,  Monte  Carlo  Methods  in  Statistical  Physics,  Springer-Yerlag,  Ber¬ 
lin.  1979 

10  D  W.  Wood,  in  Water,  A  Comprehensive  Treatise,  Vbl  6.  Recent  Advances, 
ed  F  Franks.  Plenum.  New  York.  1979,  p  279. 

II.  D  Ceperely  and  J  Tully.  Proceedings  of  the  Workshop  on  Stochastic  Molecu¬ 
lar  Dynamics.  National  Resource  for  Computation  in  Chemistry.  Berkeley, 
1979 

12  U  Levitt  and  A  Warshel.  Nature,  25a  694  (1975). 

13.  J  A  McCammon.  B  R  Gelin,  and  M  Karplus.  Nature,  067,  5B5  (1977) 

14  A.  Warshel,  in  Semi- Empirical  Methods  of  Electronic  Structure  Theory,  f\trt 
A  ed  G  Segal.  Plenum  Press.  New  York,  1977,  p  133. 

15  A  T  Hagler,  F.  Naider,  P  S  Stern,  and  R  Sharon.  J  Aon  Chem  Soc  ,  101, 
6642(1979) 

16  B  i  Alder  and  T  E  Waonwright,  J  Chem  Phys  ,  31,  459  (1969) 

17  A  Rahman,  Phys  Rev  ,  136A  405  (1964). 

18  L  Verlet,  Phys  Rev  ,  159,  9B  (1967) 

19  P  Schofield,  Cbm  Phys  Comm  ,  5,  17  (1973). 

20  J.  A  McCammon  and  M  Karplus,  Arm  Rev  phys  Chem  .,  31,  29  (1900) 

21  N  S  Ostlund,  In  Report  of  a  Minicomputer  Workshop,  National  Resource  for 
Computation  Id  Chemistry,  Berkeley,  1978. 

22  N.  S  Ostlund,  Attached  Scientific  processors  for  Chemical  Compulations:  A 
Report  to  the  Chemistry  Community,  National  Resource  for  Computation  in 
Chemistry,  Lawrence  Berkeley  Laboratory  (LBL-10409),  Berkeley,  1980,  also 
available  from  the  National  Technical  Information  Service 


£  £  %  £  S  8  $  g  8  IS  £88  88  8888 


-  22  - 


W  R  Wittmayer,  Computer  Design.  93  (March  1B?B) 

R  S  Bucy  and  K  D.  Senne,  Comp  Math  Applications.  6  317(1900) 

W  J  Karp! us  and  D  Cohen,  Computer.  11  (Sept  1931) 

D  Berg  mark,  in  Proceedings  of  the  1978  ARRAY  Conference.  Floating  Point 
Systems,  Inc.,  Portland,  1978 

Array  Processor  Fortran.  Floating  Point  Systems,  Inc  ,  Portland.  1976 

Toast  Fortran  Development  System.  System  Software  Factors,  Reading- 
Berkshire,  England,  1981 

A  E  Charlesworth,  Computer.  IB  (Sept  1961) 

H  C  Andersen,  private  communication 

C  Pottle,  M  S  Pottle,  R  W  Tuttle,  R  J  Kinch.  and  H  A  Scheraga,  J  Comp 
Chem  .  1,  46  (1980) 

D  C  Rapaport  and  H  A  Scheraga,  Chem  Ffiys  Lett  ,  78  491  (19Bl). 

R  H  Kincaid  and  H  A  Scheraga,  J  Ffiys  Chem  ,  86,  833  (1982) 

R  H  Kincaid  and  H  A  Scheraga  J  Ffiys  Chem  .  65  638(1962) 

W  C  Swope,  H  C  Andersen,  P  H  Berens.  and  K  R  Wilson.  J  Chem  Ffiys  . 

76,  637  (1982) 

B  J  Berne,  private  communication 
R  Dammkoehler,  private  communication 
P  Alexander,  Fndust  Res  Devel  ,111  (May  19B0) 

If  Furey,  Jr  .  B  C  Wang,  and  M  Sax,  J  AppA  Crystallogr  .  15  160  (19B2) 

J  M  Dawson.  R  W  Huff,  and  C  Wu.  AF1PS  National  Computer  Conference 
Proceedings.  *7,  395  (1978) 

P  H  Berens  and  K  R  Wilson,  in  Picosecond  Phenomena  II,  ed  R 
Hochstrasser.  W  Kaiser,  and  C.  V.  Shank,  Springer-Verlag.  Berlin,  19B0,  p 
246 

P.  H  Berens  and  K  R  Wilson,  J  Chem  Ffvys  ,  74,  4B72  (196:) 

P  H  Berens,  S  R  White,  and  K  R  Wilson.  J  Chem  Fftys  ,  75.  515  (198l) 

P  Bado,  P.  H  Berens,  and  K  R  Wilson,  in  Picosecond  lasers  and  Applica¬ 
tions.  ed  L  S  Goldberg,  Proc  Soc  Photo-Optic  Engin..  Bellingham.  WA. 
1962.  p  230 

45  P  H  Berens.  J  P  Bergsma,  and  K  R  Wilson,  "Molecular  dynamics  and 
spectra  III  Transient  Raman  spectra  for  a  chemical  reaction,"  to  be  sub¬ 
mitted 

46  P.  Bado,  P  H  Berens,  J  P.  Bergsma.  S  B  Wilson,  K  R  Wilson,  and  E  J 
Heller,  in  Picosecond  Ffienomena  III,  ed  K.  Eisenthal,  R  Hochstasser,  W. 
Kaiser,  and  A  laubereau.  Springer-Verlag,  Berlin.  1982,  in  press 

47  P.  H  Berens.  J  P.  Bergsma,  and  K  R  Wilson,  to  be  submitted 

48  P  H  Berens,  D  H  J  Mackay,  G  M  White,  and  X.  R  Wilson,  "Molecular 
dynamics,  thermodynamics,  and  quantum  corrections  for  liquid  water."  to 
be  submitted. 

48  B  W  Kernighan  and  D  M  Ritchie,  The  C Programming  Language,  Prentice- 
Hall.  Englewood  Cliffs,  New  Jersey,  1978. 

60  D  J.  Adams,  in  The  Problem  of  Long-Range  Forces  in  the  Computer  Simula¬ 
tion  of  Condensed  Media,  ed  D  Ceperely,  National  Resource  for  Computa¬ 
tion  in  Chemistry,  Berkeley.  1980,  p  13 


-  23- 


61.  H  V!  Cundy  and  A  P  Rollet,  Mathematical  Models,  Oxford  University  Press 
London,  1961 

62  T  Andrea,  W  C.  Swope,  and  H  C  Andersea  to  be  submitted 

63  B  Quentrec  and  C  Brot,  J  Comp  Fhys  ,  IS,  430  (1973) 

54  R  0  Watts,  Chem  Ffiys  ,  28,  367  (1977) 

55  E  B  Wilson,  Jr..  J  C  Decius,  and  P  C  Cross,  Molecular  Vibrations,  VcGravs- 
Hill,  New  York.  1955,  Chapter  4. 

56  S  Califano,  Vibrational  States.  Wiley.  London,  1976,  Chapter  4 

57  D  Beeman,  J  Comput  FTiys  .  20,  130  (1976) 

58  A  T  Hauler,  private  communication 

59  G  Chester,  R  Gann,  R  Gallagher,  and  A  Grimison,  in  Computer  Modeling  of 
Matter,  ed  P.  Lykos,  American  Chemical  Society.  Washington,  1980.  p  111 

60  A.  Komornicki,  private  communication 

61  D  Jacobs.  J  Prins,  and  K  Wilson,  in  Proceedings  of  the  1982  ARRAY  Confer¬ 
ence,  Floating  Point  Systems  Inc  .  Portland,  1982 

62  D  V.  Ritchie  and  K  Thompson,  Bell  Sys  Tech  J.,  57.  1905  (1978) 

63  R  Thomas  and  J  Yates,  A  User  CAcide  to  the  UND(*System.  Osborne,  Berkeley, 
1962 

64  P.  1  Richards,  Manual  of  Mathematical  Physics,  Pergamon.  London,  1959,  p 
299 

65  R  E  Williamson,  R  H  Crowell,  and  H  F  Trotter,  Calculus  of  Vector  Func¬ 
tions.  Prentice-Hall,  Englewood  Cliffs,  1968,  p  167. 


SP472-3/A1 


472: GAN: 716  ten j 
78u472*608 


TECHNICAL  REPORT  DISTRIBUTION  LIST,  GEN 


Office  of  Navel  Research 

Attn:  Code  472 

800  North  Quincy  Street 

No. 

Copies 

D.S.  Army  Research  Office 

Attn:  CRD-AA-IP 

P.0.  Box  1211 

No. 

Copies 

Arlington,  Virginia  22217 

ONR  Western  Regional  Office 

Attn:  Dr.  R.  J.  Marcus 

2 

Research  Triangle  Park,  N.C.  27709 

Naval  Ocean  Systems  Center 

Attn:  Mr.  Joe  McCartney 

1 

1030  East  Green  Street 

Pasadena,  California  91106 

ONR  Eastern  Regional  Office 

Attn:  Dr.  L.  H.  Peebles 

1 

San  Diego,  California  92152 

Naval  Weapons  Center 

Attn:  Dr.  A.  B.  Amster, 

Chemistry  Division 

1 

Building  114,  Section  D 

666  Summer  Street 

China  Lake,  California  93555 

1 

Boston,  Massachusetts  02210 

1 

Naval  Civil  Engineering  Laboratory 
Attn:  Dr.  R.  W.  Drisko 

Director,  Naval  Research  Laboratory 
Attn:  Code  6100 

Port  Hueneme,  California  93401 

1 

Washington,  D.C.  20390 

1 

Department  of  Physics  &  Chemistry 
Naval  Postgraduate  School 

The  Assistant  Secretary 
of  the  Navy  (RE&S) 

Department  of  the  Navy 

Room  4E736,  Pentagon 

Monterey,  California  93940 

Scientific  Advisor 

Commandant  of  the  Marine  Corps 

1 

Washington,  D.C.  20350 

Commander,  Naval  Air  Systems  Command 
Attn:  Code  310C  (H.  Rosenwasser) 
Department  of  the  Navy 

1 

(Code  RD-1) 

Washington,  D.C.  20380 

Naval  Ship  Research  and  Development 
Center 

1 

Washington,  D.C.  20360 

1 

Attn:  Dr.  G.  Bosmajlan,  Applied 
Chemistry  Division 

Defense  Technical  Information  Center 
Building  5,  Cameron  Station 

Annapolis,  Maryland  21401 

1 

Alexandria,  Virginia  22314 

Dr.  Fred  Saalfeld 

12 

Naval  Ocean  Systems  Center 

Attn:  Dr.  S.  Yamamoto,  Marine 
Sciences  Division 

Chemistry  Division,  Code  6100 

Naval  Research  Laboratory 

San  Diego,  California  91232 

I 

Washington,  D.C.  20375 

1 

Mr.  John  Boyle 

Materials  Branch 

Naval  Ship  Engineering  Center 
Philadelphia,  Pennsylvania  19112 

1 

SP472-3/A3 


i 


472:GAN:716:enj 
78u472-608 

TECHNICAL  EE PORT  DISTRIBUTION  LIST,  GEN 
No. 

Copies 

Mr.  James  Kelley 
DTNSRDC  Code  2803 

Annapolis,  Maryland  21402  1 


472 :CAN :?}  6 : tan 
78u472-608 


TECHNICAL  REPORT  DISTRIBUTION  LIST,  Q51A 

Ho. 

Copi cs 

Dr.  H.  Rauhut 
American  Cyanamid  Company 
Chemical  Research  Division 
Bound  Brook,  New  Jersey  08805  1 

Dr.  J.  I.  Zink 

University  of  California,  Los  Angeles 
Department  of  Chemistry 
Los  Angeles,  California  90024  ] 

Dr.  B.  Schechtman 
IBM 

San  Jose  Research  Center 
5600  Cottle  Road 

San  Jose,  California  95143  1 


Dr.  C.  A.  Heller 
Naval  Weapons  Center 
Code  »>059 

China  Lake,  California  93555  1 

Dr.  J.  R.  MacDonald 
Naval  Research  Laboratory 
Chcniitry  Division 
Code  hi  10 

Washington,  D.C.  20375  J 

Dr.  C.  B.  Schuster 
University  of  Illinois 
Chemi*  try  Department 

Urban..,  Illinois  61801  1 

Dr.  E.  M.  Eyring 
University  of  Utah 
Department  of  Chemistry 

Saif  Lake  City,  Utah  1 


Dr.  John  Cooper 
Code  6130 

Naval  Research  Laboratory 

Washington,  D.C.  20375  1 


if ornia , 


$4 


£bemistry 
f ornia 


90024 


No. 

Copies 


Dr.  M.  W.  Windsor 
Washirgton  State  University 
Department  of  Chemistry 
Pullnun,  Washington  99163 

Dr.  C.  R.  Bernstein 
Colori.do  State  University 
Department  of  Chemistry 
Fort  Collins,  Colorado  80521 


Dr.  A.  Adamson 

Univeisity  of  Southern  California 

Department  of  Chemistry 

Los  Argeles,  California  90007  1 

Dr.  .  S.  Wrighton 
Kj-.s  i  :usetts  Institute  of 
U«.:-  .ology 

Department  of  Chemistry 

Cambrid.  .-,  ‘lassarh  isetts  02139  1 


SP472-3/A27 


472:GAN;-;6:ddc 

7PU472-6C8 


TECHNICAL  REPORT  DISTRIBUTION  LIST ,  051? 
No. 

CoTTes 


Professor  K.  Vilipd' 

De£fc<rTBent  of  Cheaistrv , 

Univers  ijyhoiCalifornia , 

Sen  Mego 

Le  J^lla ,  California  92093 

Professor  C.  A.  Anaell 
Deoartaent  of  Chemistrv 
Purdue  University 
West  Lafayette,  indiene  47907 

Professor  P.  Meiier 
Deoertaent  of  Physics 
Catholic  University  of  America 
Washington,  D.C.  20064 

Dr.  S.  Greer 
Chemistry  Department 
University  of  waryland 
Colleee  Perk,  Maryland  20742 

Professor  ?.  Delahay 
New  York  University 
100  Washington  Souere  East 
vew  York,  New  YorV  10003 

Dr.  T.  Ashworth 
Deoertaent  of  Physics 
South  Dakota  School  of 
Mines  A  Technology 
Rapid  City,  South  Dakota  57“C1 

Dr.  G.  Gross 
New  Mexico  Institute  of 
Mining  S  Teehnoloey 
Socorro,  New  Mexico  57P.01 


Dr.  B.  Vonnegut 
State  University  of  New  York 
Earth  Sciences  Building 
1400  Washington  Avenue 
Albany,  New  York  12203 

Dr.  Rank  Loos 

Laguna  Research  Laboratory 

21421  Stans  Lane 

Laguna  Beach,  California  92651 

Dr.  John  Latham 
University  of  Manchester 
Institute  of  Science  &  Teehnoloey 
P.0.  Box  88 

Manchester,  England  M601QP 


Dr.  J.  Kassner 

Soace  Science  Research  Center 
University  of  Missouri  -  Rolla 
Rolla,  Missouri  65401 

Dr.  J.  Telford 
Universitv  of  Nevada  System 
Desert  Research  Institute 
Lab  of  Atmospheric  Physies 
Reno,  Nevada  39507 


SP472-3/B1 


472 :GAN: 716 : lab 
78u472-608 


TECHNICAL  REPORT  DISTRIBUTION  LIST,  Q51C 


No. 

Copies 

Dr.  N.  B.  Denton 
Department  of  Chemistry 
University  of  Arizona 

Tucson,  Arizona  85721  1 

Dr.  R.  A.  Osteryoung 
Department  of  Chemistry 
State  University  of  New  York 
at  Buffalo 

Buffalo,  New  York  14214  1 

Dr.  B.  R.  Kowalski 

Department  of  Chemistry 

University  of  Washington 

Seattle,  Washington  98105  1 

Dr.  S.  P.  Perone 
Department  of  Chemistry 
Purdue  University 

Lafayette,  Indiana  47907  1 

Dr.  D.  L.  Venezky 
Naval  Research  Laboratory 
Code  6130 

Washington,  D.C.  20375  1 

Dr.  H.  Freiser 
Department  of  Chemistry 
University  of  Arizona 
Tuscon,  Arizona  85721 

Dr.  Fred  Saalfeld 
Naval  Research  Laboratory 
Code  6110 

Washington,  D.C.  20375  1 

Dr.  H.  Chernoff 
Department  of  Mathematics 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02139  1 

Dr.  KT'W^itfon 

Depaprwnt  of~CHlMstry 

University  of  California,  San  Diego 

La  Jolla,  California  1 

Dr.  A.  Zirino 
Naval  Undersea  Center 

San  Dlago,  California  92132  1 


No. 

Copies 

Dr.  John  Duff in 

United  States  Naval  Postgraduate 
School 

Monterey,  California  93940  1 

Dr.  G.  M.  Hief t je 
Department  of  Chemistry 
Indiana  University 

Bloomington,  Indiana  47401  1 

Dr.  Victor  L.  Rehn 
Naval  Weapons  Center 
Code  3813 

China  Lake,  California  93555  1 

Dr.  Christie  G.  Enke 

Michigan  State  University 

Department  of  Chemistry 

East  Lansing,  Michigan  48824  1 

Dr.  Kent  Eisentraut,  MBT 

Air  Force  Materials  Laboratory 

Wright-Patterson  AFB,  Ohio  45433  1 

Walter  G.  Cox,  Code  3632 
Naval  Underwater  Systems  Center 


Building  148 

Newport,  Rhode  Island  02840  1 

Professor  Isiah  M.  Warner 

Texas  A&M  University 

Department  of  Chemistry 

College  Station,  Texas  77840  1 

Professor  George  H.  Morrison 

Cornell  University 

Department  of  Chendsty 

Ithaca,  New  York  14853  1 

Professor  J.  Janata 
Department  of  Bioengineering 
University  of  Utah 

Salt  Lake  City,  Utah  84112  1 

Dr.  Carl  Heller 
Naval  Weapons  Center 

Chins  Lske,  California  93555  1 


_  -‘iwrtllWn  i  •  - 


j 


1 


472:GAN:716:lab 

78u472-608 


TECHNICAL  REPORT  DISTRIBUTION  LIST,  Q51C 
No. 

Copies 


Dr.  L.  Jervis 
Code  A 100 

Naval  Research  Laboratory 

Washington,  D.C.  20375  1 


2 


